library(ggplot2)
library(data.table)
library(rjson)
library(groHMM)
library(dplyr)
library(cowplot)
library(pbmcapply)
library(default)
library(gganimate)
library(weights)

rm(list = ls())
source('meta.R')
source('events.R')
source('queue.R')
source('movie.R')

Inputs

inputdir <- '../data/C1D.matrix.gz' # input file is a matrix generated by computeMatrix of deepTools
configdir <- '../data/config.txt' # a configuration file for parameters
spikeinList <- c(19163158, 17625265, 17986917, 18285376, 20451950, 24631397) # length of the spikein list should be the same as the sample number

Configuration

config <- read.delim(configdir, header = F, stringsAsFactors = F, comment.char = "#")
for(i in 1:nrow(config)){
  assign(config[i, 1], config[i, 2])
}

inputMeta <- combineMeta(inputdir, spikeinList, outputdir)
metaGMerge <- inputMeta[[1]]
metaGPara <- inputMeta[[2]]
binsize <- inputMeta[[3]]

maxl <- max(metaGMerge$x) * binsize
maxPeak <- metaGPara[1, "peak"] * binsize
maxEnd <- metaGPara[1, "end"] * binsize

# reconstruct functions
default(fvx) <- list(phia = phia, phib = phib, phic = phic, keln = keln, maxPeak = maxPeak)
default(pvev) <- list(alpha = alpha, lambda = lambda, miu = (maxPeak+maxEnd)/2, sigma = (maxEnd-maxPeak)/3)
default(upho) <- list(theta = theta, gamma = gamma, taoa = taoa, taob = taob, u0a = u0a, u0b = u0b, u0c = u0c)
default(phoTrans) <- list(rho = rho)

Preload trajectory

# generate a cheat list of trajectory for saving time :)
t0 <- ifelse(tmin<=0, tmin - tstable, 0-tstable)
moveSlowest <- tmin-t0
moveTrack <- pbmclapply(seq(t0, tmax, tspan), function(ti){
  fpos(ti, (ti+moveSlowest), tspan)
}, mc.cores = ncores) %>%
  do.call(rbind, .) %>%
  as.data.frame()
rownames(moveTrack) <- as.character(round(seq(t0, tmax, tspan), 10))
colnames(moveTrack) <- as.character(round(seq(0, ceiling(moveSlowest), tspan), 10))

Recycling process after CTDP1 degradation

# recycle
ReRes <- ReSingleCell(pfree = pfree, thres = thres, moveTrack = moveTrack)

Normalize and combine track

ReTrack <- ReTrackNorm(ReRes[[1]], metaGMerge)
plotGMergeRe(ReTrack, metaGMerge)

Recycling factor

ReDf <- data.frame(t = seq(tmin, tmax, tspan), pool = ReRes[[2]], nload = ReRes[[3]], ndrop = ReRes[[4]])
ReDf$RF <- ReDf$pool / ReDf$pool[1]
ggplot(ReDf, aes(x = t, y = RF))+
  geom_line() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ylab("Recycling factor") +
  scale_x_continuous(breaks = seq(tmin, 180, by = 30))+
  theme_classic()+
  theme(axis.text = element_text(color = "black"),
        axis.ticks = element_line(color = "black"))

Create movie

It takes about 2 hours to generate a movie from -30 to 185 min with 100 cells using 15 cores.

viniSto = vini * 10
pfreeSto = viniSto * pfree
thresSto = viniSto * thres

print('Generating cells ...')
plotdfCell <- pbmcapply::pbmclapply(1:ncell, function(i){
  singleCell(aseed = i*55555)
}, mc.cores=ncores)

combineCell <- lapply(plotdfCell, function(x){x[[1]]}) %>%
  rbindlist() %>%
  split(., .$t) %>%
  mclapply(., countPol2)

combineTrack <- do.call(rbind, combineCell) %>% as.data.frame()
combineTrack$t <- as.numeric(rownames(combineTrack))
plotTrack <- reshape2::melt(combineTrack, id.vars = "t", variable.name = "pos", value.name = "y")
plotTrack$pos <- as.numeric(levels(plotTrack$pos))[plotTrack$pos]
plotTrack$t <- as.numeric(plotTrack$t)

print('Generating mean effective pool size ...')
plotRe <- lapply(plotdfCell, function(x){x[[2]]}) %>%
  rbindlist() %>%
  group_by(t) %>%
  summarise(pool = mean(pool), loadn=mean(loadn), dropn=mean(dropn), vini=mean(vini)) %>%
  as.data.frame()
plotRe$poolNorm <- plotRe$pool / pfreeSto
plotRe$pos <- 45
plotRe$y <- 2

print('Combine movement and track ...')
plotdf <- plotdfCell[[1]][[1]] # pick a cell for visulization
## [1] "Generating cells ..."
## [1] "Generating mean effective pool size ..."
## [1] "Combine movement and track ..."
anim <- ggplot(plotdf, aes(x = pos, y = y)) +
  geom_point(size = 1.5, aes(color = upho)) +
  geom_point(data = plotRe, aes(size = poolNorm), fill = "white", shape = 21, show.legend = F) +
  annotate("text", x = 45, y = 3, label = "effective Pol II pool") +
  geom_line(data = plotTrack) +
  labs(color = "pSer2 per Pol II") +
  scale_color_viridis_c(option = "magma", end = 0.9, limits=c(1,2.25)) +
  scale_size(range = c(20*min(plotRe$poolNorm), 20)  )+
  scale_x_continuous(breaks = c(0, 35, maxl), limits = c(0,maxl), labels = c("TSS","TES","TES+15kb")) +
  scale_y_continuous(limits = c(-1, 4)) +
  theme_classic() +
  theme(axis.text = element_text(color = "black"),
        axis.line = element_line(color = "black")) +
  transition_time(t) +
  labs(title = "time: {round(frame_time, 1)} min", x = NULL, y = "pSer2 state signals") +
  ease_aes('linear')
animate(anim, nframes = length(unique(plotdf$t)), height = 3, width = 8, units = "in", res = 300)
anim_save("../results/model.movie.gif")